# pests.r
#
# A PESTS source for R  (and maybe S-plus)
#
# PESTS = Poisson Estimators for State-space Time Series
# VERSION 1.1.5
#
# Requires R version >= 2.9.0
#
#
# Patrick Brandt
# pbrandt@utdallas.edu
# School of Economic, Political, and Policy Sciences
# The University of Texas, Dallas
# 800 W. Campbell Rd., GR 31
# Richardson, Texas 75080
#
# This file contains no warranty or guarantee of performance.
# Use at your own risk for non-commercial use
#
# Contents based on
#
# Brandt et al. 2000. "Dynamic Modelling for Persistent Event
# Count Time Series." American Journal of Political Science. 44(4):
# 823-843.
#
# Brandt and Williams. 2001. "A Linear Poisson Autoregressive
# Model: The PAR(p)"  Political Analysis 9(2).
#
# IF YOU USE THESE FUNCTIONS, CITE THE ABOVE REFERENCES.
#
# These functions are translated with modifications from the original
# GAUSS code.
# Copyright 1998-2005, Patrick T. Brandt and John T. Williams.
#
# R version of the code Copyright 2001-2006, Patrick T. Brandt.
#
# Revision history
#
# 20061120
#
# 20061025 Added descriptions of the functions for the PAR(p)
#          multiplier functions in this file
#
# 20050920 Revised the Pewma() function to estimate models WITHOUT any
#          covariates or a model with ONLY an intercept.
#
# 20041021 Added options to parp() and Pewma() for user supplied
#          starting values (thanks to Dan Wood for requesting this).
#
# 20040419 Revised the function calls for PEWMA and PAR(p) to use the
#          formula interface.
#          Implemented a class space for the models.
#          Implemented a new output format for the results.
#          Include a default test of omega=1 for the PEWMA and rho_i =
#          0 for PAR(p).
#          Added print functions for the class of the models.
# 20090907 Corrected some of the logic for how user input initial
#          starting values for the pewma() function are handled.
#---------------------------------------------------------------------#


####################################################################
# EWMA function for use in the PEWMA model
# Computes a simple EWMA
# Inputs:
# x = vector
# omega = weight parameter in [0.1]

ewma<-function(x,omega)
{ ind<-omega^(rev(seq(1:length(x))))
  wm<-weighted.mean(x,ind)
  wm
}

####################################################################
# Purpose: Implements a Wald Test using the covariance matrix and
# a restriction matrix for general OLS / MLS regressions
#
# Format: wald(coef,restrict,covar,t);
#
# INPUTS: coef = Kx1, vector, coefficent vector
#     restrict = RxK, matrix, restriction matrix
#     covar = KxK, matrix, covariance matrix (from estimation)
#          that is symmetric and PD
#     t = scalar, sample size for df calculation
#
# OUTPUTS: wald_val = Wald statistic under the null
#          wp = p-value for X^2 with R degrees of freedom

wald<-function(coef,restrict,covar,t)
{ # Trap non-conformable restrictions */
  if(ncol(restrict)!=nrow(coef))
  {
    stop("WALD TEST ERROR: Restriction and coefficient matrices are non-conformable")
  }
  # Compute Wald statistic and p-value
  waldval<-t(coef)%*%t(restrict)%*%solve(restrict%*%covar%*%t(restrict))%*%restrict%*%coef
  # Small sample correction
  waldval <- ((t-nrow(coef)+nrow(restrict))/t)*waldval
  wp<-1-pchisq(waldval,nrow(restrict));
  
  list(Wald.stat=waldval,Wald.p=wp);
}

####################################################################
# PEWMA code
# There are 4 functions:
# 1) DGP function
# 2) PEWMA Kalman filter
# 3) PEWMA model log-likelihood function
# 4) PEWMA wrapper that estimates a PEWMA regression model using
#    optim()
####################################################################

####################################################################
# pewmadgp -- PEWMA data generation process
# n = number of observations in generated series
# a0 = initial value of a0, shape parameter for gamma prior
# b0 = initial value of b0, scale parameter for gamma prior
# w = omega, EWMA weight parameter
# x = matrix of regressors, NOT including an intercept, unless you
#     want a model with a growth term.
# d = delta, vector of regression coefficients for x.
#

pewmadgp<-function(n,a0,b0,w,x,d=c(0))
{ y<-matrix(0,nrow=n); mu<-matrix(0,nrow=n);
  d<-matrix(d,ncol=1)
  r<-matrix(0, nrow=n, ncol=1);
  for(t in 1:n)
  {   if(t==1)
  { r[t]<-digamma(a0) - digamma(w*a0)
    m<-(a0*w*exp(r[t]))/(b0*w*exp(-x[t,]%*%d))
    y[t]<-rpois(1,m)
    # update to compute posterior
    a<-a0*w + y[t]
    b<-b0*w + exp(x[t,]%*%d + r[t]);
    mu[t]<-y[t]
  }
  else
  { r[t]<-digamma(a) - digamma(w*a)
    aa<-a*w
    bb <- b*w*exp(-x[t,]%*%d - r[t])
    m<-(exp(x[t,]%*%d + r[t])*(ewma(y[1:t-1],w))/(ewma((exp(x[1:t-1,]%*%d + r[1:t-1])),w)))
    y[t]<-rpois(1,m)
    rm(m)
    a<-aa + y[t]
    b<-w*b + exp(x[t,]%*%d + r[t]);
    mu[t]<-exp(x[t,]%*%d + r[t])*a/b;
  }
  }
  y<-ts(y,start=1)
  y
}

####################################################################
# pewmafilter -- Extended Kalman filter for PEWMA model.  This is used
# to compute the filter and the parameters for the log-likelihood
# function, pewmallf().
# THIS FUNCTION IS NOT CALLED BY THE USER -- IT IS USED IN THE MAIN
# Pewma() ESTIMATION FUNCTION

pewmafilter<-function(y,x,w,d)
{  # Initialize some matrices we need to store the filter
  a<-matrix(0,nrow=length(y)); aa<-a; aaa<-a;
  b<-a; bb<-a; bbb<-a;
  mu<-a; r<-a;
  # Loop over the filter calculations
  for(t in 1:length(y))
  { if(t==1)
  { aaa[t]<-y[t]; bbb[t]<-exp(x[t,]%*%d + r[t,1]);
    mu[t]<-mean(y);
  }
  else
  { a[t]<-aaa[t-1];  b[t,1]<-bbb[t-1];
    r[t]<-digamma(a[t]) - digamma(w*a[t]);
    aa[t]<-w*a[t,1];
    bb[t]<-w*b[t,1]*exp(-x[t,]%*%d)*exp(-r[t]);
    aaa[t]<-w*a[t,1] + y[t];
    bbb[t]<-w*b[t,1] + exp(x[t,]%*%d + r[t]);
    mu[t]<-exp(x[t,]%*%d + r[t])*aaa[t]/bbb[t];
  }
  }
  # Clean up and store the filter components in a ts object
  dta<-matrix(cbind(a,aa,aaa,b,bb,bbb,mu,r),ncol=8)
  dta<-ts(dta, start=1, names=c("a","aa","aaa","b","bb","bbb","mu","r"))
  dta
}

####################################################################
# pewmallf -- PEWMA log-likelihood function.
# Returns the value of the LLF for the given data and parameters.

pewmallf<-function(p,y=y,x=matrix(0,nrow=length(y)))
{ # initialize the parameters
  if(length(p)==1)
  { w<-p[1];  d<-as.matrix(0) }
  else
  { w<-p[1]; d<-as.matrix(p[2:length(p)]) }
  y<-as.ts(y); x<-as.ts(x);
  # Call the filter
  pf<-pewmafilter(y,x,w,d);
  # Extract vectors of parameters from filter output
  pf<-pf[2:nrow(pf),];
  a<-pf[,"a"]; b<-pf[,"b"]; r<-pf[,"r"];
  x<-x[2:nrow(x),];    y<-y[2:nrow(y)];
  # Compute log-likelihood function
  lp<-(lgamma((w*a) + y) - lgamma(w*a)
       - lgamma(y+1) + ((w*a) * log(w*b*exp(-x%*%d - r)))
       - ((w*a) + y)*(log(1 + (w*b*exp(-x%*%d - r)))));
  lp<-sum(lp)
  lp
}

####################################################################
# Pewma -- Fits a Poisson Exponentially Weighted Moving Average Model
#          by maximum likelihood.
#
# Syntax: Pewma(formula, omega.init=0.6, init.param=NULL, ...)
#
# Description:
#
# This function uses the standard formula interface:
# y ~ x1 + x2 + ...
#
# y = Event count time series you wish to fit.  (can be a matrix or a
# ts object)
# x = regressors not including an intercept, unless you
# want a growth term.
#
# init.param = starting values for the optimizer.  Use a glm() call to
# compute these or leave NULL and they will be computed for you.
#
# THIS IS IMPORTANT!  IF YOU DO *NOT* THINK THERE
# IS EXPONENTIAL GROWTH IN THE COUNTS USE A FORMULA OF THE FORM:
#
# y ~ -1 + x1 + x2 ....
#
# If you think the counts are a random walk without drift you want the
# above specification, which estimates a model that has a latent
# random walk with zero drift. (IF YOU DO NOT KNOW WHAT THIS IS
# SAYING, USE THE ABOVE SPECIFICATION WITH THE "MINUS 1" TERM IN THE
# FORMULA!)
#
# If there is random drift, you want an intercept or a model formula like
# this:
#
# y ~ 1 + x1 + x2 ....
#
# This model will have drift in the latent mean.
#
# Models without covariates can now be fit using the following formula
# specifications:
#
# y ~ -1  : Just the PEWMA process -- estimates only omega.
# y ~ 1   : PEWMA process with drift -- estimates omega and a trend
#           term as the intercept
#
# Other optional arguments:
#
# omega.init = 0.6 -- starting value for omega.  Can be changed, but
#                     0 < omega <= 1
# init.param = NULL --- Function uses a Poisson regression for
#                       starting values.  This can be replaced with a
#                       list, c(....) for user supplied starting
#                       values.  Must be conformable.



# Note that the R version uses the constrained optimizer for the
# parameter w.  This is because 0 < w <= 1.  The parameter can go
# slightly over the bound of 1 because this is where the Pewma reduces
# to the Poisson model.  If the bound bites for d, just alter the
# value of 1000 to something else.

Pewma<-function(formula, omega.init=0.6, init.param=NULL, ...)
{ # Set up the model
  
  model <- model.frame(formula)
  
  # Setup of x and y for the model depends on the number of
  # covariates and the inclusion of the intercept
  
  # Get y and x
  y <- model.extract(model, "response")
  x <- as.matrix(model.matrix(formula))
  xnames <- dimnames(x)[[2]]
  k <- ncol(x)
  
  # Get starting values if needed
  if(is.null(init.param)==TRUE & k>0)
  {
    init.param <- glm(formula, family=poisson())$coefficients
  }
  
  # Now build the correct bounds for the optimizer
  
  bound <- 1000
  if(k>0)
  { lower.bound <- c(0,rep(-bound,k))
    upper.bound <- c(1.1,rep(bound,k))
  } else {
    lower.bound <- c(0)
    upper.bound <- c(1.1)
  }
  
  # Fit the model
  if(k>0)
  {
    fitted.param<-optim(c(omega.init,init.param), pewmallf, method = "BFGS",
                        #                            lower=lower.bound, upper=upper.bound,
                        hessian=TRUE,
                        control=list(trace=2, REPORT=5, fnscale=-(length(y))),
                        y=as.matrix(y),x=as.matrix(x))
  } else {
    fitted.param<-optim(omega.init, pewmallf, method = "L-BFGS-B",
                        lower=lower.bound, upper=upper.bound,
                        hessian=TRUE,
                        control=list(trace=2, REPORT=5, fnscale=-(length(y))),
                        y=as.matrix(y))
  }
  
  # Retrieve the parameters and the llf value
  param<-fitted.param$par
  llf <- fitted.param$value
  
  # Estimate the standard errors and compute the z scores
  if(k>0)
  { covar <- -solve(fitted.param$hessian)
    se<-sqrt(diag(covar))
    z<-param/se
  } else {
    covar <- -1/(fitted.param$hessian)
    se<-sqrt((covar))
    z<-param/se
  }
  coefs <- t(matrix(rbind(param,se,z),nrow=3))
  colnames(coefs) <- c("Parameters","Std. Errors","Z-score")
  
  if(k>0)
  {
    rownames(coefs) <- c("Omega", colnames(x))
  } else {
    rownames(coefs) <- c("Omega")
  }
  
  cat(" ","\n")
  cat("--------------------------------------","\n")
  cat("PEWMA regression output","\n")
  cat("--------------------------------------","\n")
  printCoefmat(coefs)
  cat("--------------------------------------","\n")
  cat("Log-likelihood value  : ", llf, "\n")
  aic <- -2*llf+2*k
  cat("AIC                   : ", aic, "\n")
  dof <- length(y)-k
  cat("Degrees of Freedom    : ", dof, "\n")
  cat("--------------------------------------","\n")
  cat("Test for reduction to Poisson\n")
  cat("--------------------------------------","\n")
  wald.test <- ((param[1]-1)/se[1])^2
  wald.p <- 1-pchisq(wald.test,1)
  cat("Wald test for omega=1 : ", wald.test, "\n")
  cat("Wald p-value          : ", wald.p, "\n")
  cat("--------------------------------------","\n\n")
  
  fit <- list(coefs=coefs,
              optim.param=fitted.param,
              param=param,
              covar=covar,
              std.err=se,
              z=z,
              llf=llf,
              aic=aic,
              dof=dof,
              wald.test=wald.test,
              wald.p=wald.p,
              k=k,
              y=y,
              x=x,
              call=match.call(),
              obj.type=c("PEWMA regression output"))
  class(fit) <- c("pests.output")
  return(fit)
}

####################################################################
# PAR(p) code
# There are 4 functions:
# 1) DGP function
# 2) PAR(p) Kalman filter
# 3) PAR(p) model log-likelihood function
# 4) PAR(p) wrapper that estimates a PAR(p) regression model using
# optim()
####################################################################

####################################################################
# parpdgp -- PAR(p) data generation process
# n = number of observations in generated series
# r = AR(p) coefficients
# x = matrix of regressors, including an intercept.
# d = delta, vector of regression coefficients for x.
#

parpdgp<-function(n,r,x=matrix(0,nrow=n),d=as.matrix(0))
{ # Initialize the parameters
  y<-matrix(0,nrow=n)
  r<-as.matrix(r)
  p<-nrow(r)
  xd<-exp(as.matrix(as.matrix(x)%*%as.vector(d)))
  # Loop over DGP
  for(t in 1:n)
  { if(t<(p+1))
  { y[t]<-rpois(1,xd[t]) }
  else
  { ylag<-t(y[(t-p):(t-1)])%*%rev(r)
    y[t]<-rpois(1,(ylag + (1-sum(r))*xd[t]))
  }
  }
  return(ts(y,start=1))
}

####################################################################
# parpfilter -- Extended Kalman filter for PAR(p) model.  This is used
# to compute the filter and the parameters for the log-likelihood
# function, parpllf().

parpfilter<-function(y,x,r,d)
{ # Initialize some parameters
  y<-as.matrix(y)
  r<-as.matrix(r)
  p<-nrow(r)
  n <- length(y)
  xd<-exp(as.matrix(as.matrix(x)%*%as.vector(d)))
  m<-matrix(0,nrow=n)
  s<-m
  # Filter loop
  for(t in 1:n)
  { if(t<(p+3))
  { m[t]<-mean(y[1:(p+3)]); s[t]<-var(y[1:(p+3)]);
  }
  else
  { ylag<-t(y[(t-p):(t-1)])%*%rev(r)
    m[t]<-ylag + (1-sum(r))*xd[t]
    s[t]<-(sum(r^2)*var(y[1:t-1]) + (1-sum(r^2))*var(xd[1:t-1]))
  }
  }
  return(ts(matrix(cbind(m,s),ncol=2),start=1,names=c("m","s2")))
  
}

####################################################################
# parpllf -- PAR(p) log-likelihood function.
# Returns the value of the LLF for the given data and parameters.

parpllf<-function(p,y=y,x=matrix(0,nrow=length(y)))
{ # initialize the parameters
  k<-ncol(x)
  no.param<-length(p)
  r<-as.matrix(p[1:(no.param-k)])
  d<-as.matrix(p[(no.param-k+1):length(p)])
  y<-as.matrix(y); x<-as.matrix(x);
  # Call the filter
  pf<-parpfilter(y,x,r,d);
  # Extract vectors of parameters from filter output
  pf<-pf[(no.param-k+1):nrow(pf),];
  m<-pf[,"m"]; s2<-pf[,"s2"];
  x<-x[(no.param-k+1):nrow(x),];    y<-y[(no.param-k+1):nrow(y)];
  # Compute log-likelihood function
  lp<-(lgamma((s2*m) + y) - lgamma(s2*m) - lgamma(y+1)
       + (s2*m*log(s2)) - ((s2*m + y)*log(1+s2)))
  lp<-sum(lp)
  lp
}

####################################################################
# Parp -- this is the main function for calling the PAR(p) model
# estimation.

# Inputs
# y = Event count time series you wish to fit.  (can be a matrix or a
# ts object)
# x = matrix of regressors,including an intercept
# formula = R formula format,, typically y ~ x1 + x2 + ...

# parp.init = REQUIRED starting the AR(p) values for the optimizer.  This is
#             the p AR parameters
# init.param = NULL.  This uses the Poisson regression values for
#             starting values unless the user supplies values in the form
#             c(d_1,d_2,...,d_k), where d_1 is the constant.
#

Parp<-function(formula, p=1, parp.init=rep(0.1,p), init.param=NULL, ...)
{
  
  model <- model.frame(formula)
  y <- model.extract(model, "response")
  x <- model.matrix(formula)
  k <- ncol(x)
  xnames <- dimnames(x)[[2]]
  colnames(x) <- xnames
  
  # Do some work to get a vector of starting values.
  if(is.null(init.param))
  {
    init.param <- glm(formula, family=poisson())$coefficients
  }
  
  # Do a simplex run to get some good starting values
  start.param <- optim(c(parp.init,init.param), parpllf,
                       method="Nelder-Mead", hessian=FALSE,
                       control=list(maxit=800, reltol=0.0001,
                                    trace=1, fnscale=-(length(y))),
                       y=y,x=x,...)
  
  # Estimate the PAR(p) model
  fitted.param<-optim(start.param$par, parpllf, method = "BFGS", hessian=TRUE,
                      control=list(maxit=1000,
                                   trace=2, fnscale=-(length(y))),
                      y=y,x=x,...)
  
  # Retreive the parameters and llf values
  param<-fitted.param$par
  llf <- fitted.param$value
  
  # Estimate the standard errors and compute z scores
  covar <- -solve(fitted.param$hessian)
  se<-sqrt(diag(covar))
  z<-param/se
  
  k <- ncol(x)
  k.plus.p <- length(se)
  num.p <- k.plus.p - k
  fitted.par.filter <- parpfilter(y,x,param[1:num.p],param[(num.p+1):k.plus.p])
  
  residuals <- y - fitted.par.filter[,"m"]
  # Print the results
  coefs <- t(matrix(rbind(param,se,z),nrow=3));
  rownames(coefs) <- c(rep("rho",num.p),colnames(x))
  colnames(coefs) <- c("Parameters","Std. Errors","Z-score")
  ##     cat(" ","\n")
  ##     cat("PAR(p) regression output","\n")
  ##     cat("--------------------------------------------","\n")
  ##     printCoefmat(coefs)
  ##     cat("--------------------------------------------","\n")
  ##     cat("Log-likelihood value  : ", llf, "\n")
  aic <- (-2*llf)+(2*(length(param)-1))
  ##    cat("AIC                   : ", aic, "\n")
  dof <- length(y) - k.plus.p
  ##    cat("Degrees of Freedom    : ", dof, "\n")
  ##    cat("--------------------------------------------","\n")
  R <- diag(1, nrow=p, ncol=k.plus.p)
  wald.test <- wald(matrix(param, ncol=1), R, covar, length(y))
  wald.p <- wald.test$Wald.p
  wald.test <- wald.test$Wald.stat
  ##     cat("Test for reduction to Poisson regression\n")
  ##     cat("--------------------------------------------","\n")
  ##     cat("Wald test for rho_i=0 : ", wald.test, "\n")
  ##     cat("Wald p-value          : ", wald.p, "\n")
  ##     cat("--------------------------------------------","\n\n")
  fit <- list(coefs=coefs,
              optim.param=fitted.param,
              param=param,
              covar=covar,
              std.err=se,
              z=z,
              residuals=residuals,
              llf=llf,
              aic=aic,
              dof=dof,
              wald.test=wald.test,
              wald.p=wald.p,
              k=k,
              p=num.p,
              y=y,
              x=x,
              call=match.call(),
              obj.type=c("PAR(p) regression output"))
  class(fit) <- c("pests.output")
  return(fit)
}


# Print function for the PEWMA and PAR(p) output.

"print.pests.output" <- function(x)
{ cat(" ","\n")
  cat(x$obj.type, "\n")
  cat("\nCall: ", deparse(x$call), "\n\n")
  cat("--------------------------------------------","\n")
  printCoefmat(x$coefs)
  cat("--------------------------------------------","\n")
  cat("Log-likelihood value  : ", x$llf, "\n")
  aic <- (-2*x$llf)+(2*(length(x$param)-1))
  cat("AIC                   : ", x$aic, "\n")
  cat("Degrees of Freedom    : ", x$dof, "\n")
  cat("--------------------------------------------","\n")
  cat("Test for reduction to Poisson\n")
  cat("--------------------------------------------","\n")
  cat("Wald test statistic   : ", x$wald.test, "\n")
  cat("Wald p-value          : ", x$wald.p, "\n")
  cat("--------------------------------------------","\n")
  
}

# Holder function for the summary method.
"summary.pests.output" <- print.pests.output

####################################################################
# ACF function standardized for count data...
# Feed in your counts and you can generate some ACFs!

# Computes an ACF for counts by standardizing them.  This is a valid
# method for searching for autocorrelation in counts.  See Cameron and
# Trivedi (Regression Analysis of Count Data, CUP, 1998: 228).

counts.acf <- function(y, lag.max=(length(y)/4),...)
{ z <- y-mean(y)
  ct.acf <- acf(z, lag.max=lag.max,plot=F,...)
  return(ct.acf)
}


# Box-Ljung Q-text for count data.  Tests for serial correlation ag
# "lag" value.
counts.Q <- function(y, lag=1, ...)
{ z <- y-mean(y)
  test <- Box.test(z,lag=lag,type=c("Ljung-Box"))
  return(test)
}


# Multiplers functions
# These work, but at present are undocumented.


# parp.multipliers() Computes impact multipliers for the PAR(p) model
# output.
#
# Inputs: parp.obj : output of a Parp() model call.  See Parp() above
#                    for details.
#         xvec : vector of values at which to compute the
#                multipliers.  This should be conformable with the number of
#                exogenous variable coefficients in the model (including the
#                constant).  TYPICALLY, this will be a vector of zeros
#                and ones, where the 1's are for the one unit change in a variable,
#                and the zeros leave a variable unchanged.  Details are documents in
#                the Brandt and Williams (2001) paper in the section on
#                interpretation of the PAR(p) model.  Note that code
#                below exactly matches equations 10 and 11

parp.multipliers <- function(parp.obj, xvec=matrix(rep(1,parp.obj$k),1,parp.obj$k))
{ p <- parp.obj$p
  k <- parp.obj$k
  ar.coef <- parp.obj$param[1:p]
  x.coef <- matrix(parp.obj$param[(p+1):(p+k)], ncol=1)
  sum.ar.coef <- sum(ar.coef)
  impact.multiplier <- ((1-sum.ar.coef)*exp(xvec%*%x.coef)*x.coef[,1])
  longrun.multiplier <- impact.multiplier/(1-sum.ar.coef)
  
  # Format
  impacts <- rbind(impact.multiplier, longrun.multiplier)
  colnames(impacts) <- colnames(parp.obj$x)
  rownames(impacts) <- c("Short Run", "Long Run")
  output <- list(impacts=impacts, call=parp.obj$call, changes=xvec)
  class(output) <- c("pestsmultipliers")
  return(output)
}

# print.pestsmultipliers: printing function for parp.multipliers

"print.pestsmultipliers" <- function(x)
{ print(x$impacts)
}

########################################################################
# montecarlo.parp.multipliers()
# Draws a parametric Monte Carlo sample from the PAR(p) multipliers.
# Inputs: first two are the same as in parp.multipliers.
#         n : number of samples to draw.

# Output: array of impact and long run multiplers.  Inner array
# indices are n x k where n is the number of samples and k is the
# number of variables over which the impacts are computed (i.e.,
# number of exogenous variables).  Outermost array indices separate on
# the short versus total impact.  The first array is the impact, the
# second the long run.

montecarlo.parp.multipliers <- function(parp.obj, xvec, n)
{ p <- parp.obj$p
  k <- parp.obj$k
  chol.covar <- t(chol(-solve(parp.obj$optim.param$hessian)))
  coef <- parp.obj$param
  mean.y <- mean(parp.obj$y)
  multipliers <- array(0, dim=c(n,k,2))
  
  for (i in 1:n)
  { coef.i <- coef + chol.covar%*%matrix(rnorm(length(coef)),nrow=(p+k))
    ar.coef <- coef.i[1:p]
    x.coef <- coef.i[(p+1):(p+k)]
    sum.ar.coef <- sum(ar.coef)
    # impact
    multipliers[i,,1] <- (1-sum.ar.coef)*exp(xvec%*%x.coef)*x.coef
    # long run
    multipliers[i,,2] <- multipliers[i,,1]/(1-sum.ar.coef)
  }
  return(multipliers)
}

#######################################################################
# Same as the montecarlo.parp.multipliers.  Just computes them for a
# Poisson model estimated by glm().  Typically used for comparison
# with the above.

glm.poisson.multipliers <- function(psn.glm,xvec,n)
{ chol.covar <- chol(vcov(psn.glm))
  coef <- coefficients(psn.glm)
  xvec <- matrix(xvec,nrow=1)
  impact <- matrix(0, nrow=n, ncol=length(coef))
  for (i in 1:n)
  { coef.i <- coef + chol.covar%*%matrix(rnorm(length(coef)),nrow=(length(coef)))
    impact[i,] <- t(rep(exp(xvec%*%coef.i),length(coef.i))*coef.i)
  }
  return(impact)
}


#######################################################################
# parp.impulse: Function to compute impulse responses for shocks in X
# variables for PAR(p) model.  Traces these out over a horizon of
# 1...n

# Inputs: parp.obj : output of a call to Parp() for a model
#         x.impulses : size of the shocks in the X variables.  Should
#                      probably only use this to shock one at a time!
#                      This will typically be a vector of 0-1 values
#                      where the 1's enter for the variable whose
#                       shock you want to trace.
# Output: impulse response vector.  Plot this over time to trace the
# response.

parp.impulse <- function(parp.obj, x.impulses, n)
{ p <- parp.obj$p
  k <- parp.obj$k
  ar.coef <- parp.obj$param[1:p]
  x.coef <- parp.obj$param[(p+1):(p+k)]
  mean.y <- mean(parp.obj$y)
  sum.ar.coef <- sum(ar.coef)
  impulse <- matrix(0, nrow=(n+p))
  for (i in (1+p):n)
  {
    impulse[i] <- (impulse[(i-1)])%*%ar.coef[1] + (1-sum.ar.coef)*exp(x.impulses[i,]%*%x.coef)
    if (p>1)
    { for (j in 2:p)
    { impulse[i] <- impulse[i] + impulse[i-j] %*% ar.coef[j]
    }
    }
  }
  return(impulse[2:n])
}
